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Galerkin’s  method  is  applied  to  reduce  a  high-order  PDE  system  to  a  lower-order  ODE 
system  using  POD  eigen-basis  to  develop  a  reduced  order  model  (ROM)  describing  the 
combustion  response  to  acoustic  excitations.  A  one-dimensional  reaction-advection  scalar 
equation  is  used  as  a  representative  equation  to  investigate  the  overall  approach.  Both  linear 
and  nonlinear  model  equations  are  employed  to  build  the  reduced  order  model.  The  influence 
of  different  discretization  methods  and  time  steps  on  the  ROM  are  evaluated.  Parametric 
explorations  are  performed  to  investigate  the  capabilities  of  the  lower-order  model  for 
predicting  combustion  response. 


Nomenclature 


PDE  = 

Partial  Differential  Equatioin(s) 

ODE  = 

Ordinary  Differential  Equation(s) 

ROM  = 

Reduced  Order  Model 

POD  = 

Proper  Orthogonal  Decomposition 

SVD  = 

Singular  Value  Decomposition 

PTE  = 

Elame  Transfer  Eunction 

EDE  = 

Elame  Describing  Eunction 

T 

temperature 

4(0  = 

POD  temporal  coefficient 

= 

POD  singular  value 

a(X)  = 

POD  temporal  mode 

(Hx)  = 

POD  spatial  mode/POD  eigen-basis 

N,  = 

number  of  POD  eigen-basis  included 

^p, total  — 

number  of  POD  eigen-basis  in  total 
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NI 


number  of  grid  points 
imaginary  unit 


j 

Superscripts 

'  =  perturbation  variables 

=  averaged/steady- state  variables 

Subscripts 

k,  n,  m=  POD  mode  # 

i  =  #  of  grid  point 

ave  =  time-averaged  value 


I.  Introduction 

Combustion  dynamics  are  a  complex  phenomena  that  results  from  the  coupling  between  the  modes  of  heat  release 
and  acoustics.  In  practical  combustor  devices  the  complexity  is  greatly  amplified  by  turbulent,  compressible  flow, 
very  high  rates  of  heat  release,  and  often  complicated  geometries  and  acoustic  boundary  conditions.  Modern 
computational  capability  offers  the  potential  for  moving  beyond  the  empirically -based  design  analysis  of  the  past,  but 
simulations  of  full  scale  dynamics  for  engineering  analysis  are  still  far  out  of  reach.  However,  high  fidelity  simulations 
of  smaller  scale  domains  can  yield  reduced  order  models  of  the  combustion  response  that  can  presumably  provide 
accurate  descriptions  of  the  linear/nonlinear  coupling  between  acoustics  and  combustion. 

Demands  for  developing  efficient  methodologies  to  establish  lower-order  engineering  models  that  can  describe 
the  linear/nonlinear  coupling  between  acoustics  and  combustion  are  increasing.  A  hierarchy  of  integrated  experimental 
and  computational  efforts  for  developing  acceptable  models  is  shown  diagrammatically  in  Fig.  1  proposed  by  Portillo 
et  al.^  The  key  idea  is  that  concurrent  experiments  and  computations  of  combustion  instability  phenomena  are  used  to 
provide  high  fidelity  outputs  from  which  response  functions  (flame  models)  can  be  derived  for  use  in  lower  fidelity 
engineering-level  testbeds.^  In  turn,  these  engineering  models  can  be  used  to  provide  quick  and  reasonable  predictions 
of  the  necessary  physics  of  interest  (e.g.,  growth  rate  and  limit  cycle)  in  actual  combustor  configurations  and,  thereby, 
aid  in  their  design  evolution.  Experiments  provide  direct/indirect  measurements  of  physics  but  they  are  limited  in  the 
amount  of  information  that  can  be  output  and  may  also  have  issues  in  the  spatial  and  temporal  resolution  of  the 
important  physics.  Therefore,  deriving  response  functions  from  experiments  can  be  restrictive.  On  the  other  hand, 
high-fidelity  physics-based  computations  such  as  large  eddy  simulations  (LES)  or  detached  eddy  simulations  (DES) 
can  provide  a  finely  resolved  description  of  the  physics  in  the  combustor.  Therefore,  they  form  the  ideal  basis  for  the 
derivation  of  reduced-order  models  (ROM),  which  is  the  subject  of  the  present  paper.  We  point  out  that  that  the 
experimental  measurements  are,  of  course,  valuable  as  a  means  of  validating  both  the  high-fidelity  as  well  as  the  newly 
developed  lower-fidelity  models. 

The  basic  idea  that  we  are  exploring  in  this  paper  involves  using  proper  orthogonal  decomposition  (POD) 
techniques  to  derive  a  reduced  equation  set.^  The  POD  is  applied  on  the  detailed  simulation  data  to  obtain  a  set  of 
spatial  basis  vectors,  which  are  then  used  within  a  Galerkin  formulation  that  reduces  the  original  partial  differential 
equations  to  an  ordinary  differential  equation.  The  Galerkin  method  has  been  widely  used  for  model  reduction. A 
simple  example  can  be  found  in  solving  the  classical  wave  equation,  where  Fourier  series  expansions  are  used  as  the 
eigen-bases.  Instead  of  solving  the  PDE  system  in  both  time  and  space,  after  reduction,  the  system  only  needs  to  be 
solved  in  time.  In  real  engineering  problems,  POD  eigen-bases  (also  known  as  POD  modes)  have  been  widely  used  in 
computational  and  experimental  studies  of  flow  fields  and  dynamics. In  addition,  the  application  of  POD  eigen- 
bases  for  reacting  flow  investigations  has  also  received  attention.  More  importantly,  POD  eigen-bases  have  been 

found  to  be  useful  in  model  reduction  of  cold  flow  problems.  Moreover,  reduced  order  models  derived  from 

incompressible  flow  simulations  have  been  applied  in  flow  control  problems  with  classical  control  theories. 
Rowley  et  al.^^  applied  the  POD-Galerkin  method  to  compressible  flow.  Application  of  reduced  order  models  to 
combustion  problems  has  also  been  performed  using  energy  and  species  equations.  Munipalli  et  al.^^  recently  devised 
approaches  for  applying  model  reduction  techniques  in  detailed  simulations  of  combustion  dynamics. 
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Fig.  1  Hierarchy  of  integrated  analysis,  computations  and  experimental  work.^ 


Specifically,  these  ideas  are  also  related  to  the  so-called  flame  transfer  function  (FTF)  or  flame  describing  function 
(FDF)  described  in  Ref  20  and  widely  used  in  turbulent  premixed  flame  studies  in  gas  turbine  engines. The  main 
difference  between  the  FTF/FDF  and  the  ROM  approaches  is  that  FTF/FDF  characterizes  the  system  response  based 
on  certain  inputs  (e.g.,  velocity  or  pressure  perturbations)  and  outputs  (e.g.,  flame  perturbations)  using  a  black-box 
treatment  (see  Ref  24).  Consequently,  the  responses  built  from  FTF/FDF  can  lose  generality  and  are  limited  to  certain 
operating  conditions,  combustor  configurations,  target  frequencies,  range  of  amplitudes,  which  means  that,  for  each 
specific  configuration  or  conditions,  the  FTF/FDF  would  need  to  be  derived,  tested  and  validated.  The  ROM  approach 
provides  a  strong  mathematical  framework  for  such  a  development  that  has  the  potential  to  be  generalized  to  broader 
conditions  and  configurations. 

The  specific  objective  of  this  paper  is  to  demonstrate  a  fundamental  procedure  for  deriving  a  reduced-order  model 
formulation  using  a  Galerkin  approach  based  on  proper  orthogonal  decomposition  (POD)  techniques  for  representing 
the  spatial  eigen-bases  of  the  reduced  solution.  A  scalar  one-dimensional  reaction-advection  model  equation  is  used 
in  these  studies  because  it  captures  many  of  the  essential  attributes  of  the  combustion  problem  within  an  efficient  test 
platform.  The  approach  is  to  obtain  a  series  of  numerical  solutions  of  the  governing  equation  by  perturbing  quantities 
of  interest  such  as  the  inlet  conditions  and  to  use  the  numerical  solutions  to  build  a  POD  eigen-basis.  The  POD-derived 
eigen-basis  is,  in  turn,  used  in  conjunction  with  a  Galerkin  procedure  to  reduce  the  governing  partial  differential 
equation  to  an  ordinary  differential  equation  which  constitutes  the  ROM.  Once  the  ROM  is  established,  it  can  be  used 
as  a  lower-order  test-bed  to  predict  detailed  results  within  certain  parametric  ranges  at  a  fraction  of  the  cost  of  the  full 
governing  equations.  Moreover,  the  ROM  has  the  further  potential  to  derive  reduced-order  response  functions,  which 
can  themselves  be  used  in  the  context  of  other  reduced-physics  models.  The  present  paper  is  focused  on  the 
development  and  validation  of  the  ROM.  Application  to  response  function  derivation  will  be  considered  in  future 
work. 

The  paper  is  organized  as  follows.  In  Section  II,  we  present  the  reaction-advection  equation  used  in  the  present 
study  as  well  as  the  Galerkin  formulation  and  POD  techniques  for  deriving  reduced  order  models  for  linear  and  non¬ 
linear  versions  of  the  equation.  In  Section  III,  we  describe  the  test  problem  that  is  used  in  the  validation  studies.  In 
Section  IV,  we  present  simulation  results  that  demonstrate  the  capabilities  of  the  reduced-order  model  to  predict 
“combustion”  response.  Specifically,  both  single -frequency  and  broadband-forcing  studies  are  employed  in  making 
these  evaluations.  In  the  final  section,  we  provide  some  concluding  remarks  and  future  directions  for  continued 
research  and  development. 


American  Institute  of  Aeronautics  and  Astronautics 


3 


II.  Formualtion 


A.  Model  reaction-advection  equation 

The  heat-addition  zone  in  combustion  problems,  the  heat-addition  zone,  where  the  instability  driving  mechanisms 
typically  originate,  is  characterized  by  a  region  with  exponential  temperature  rise  followed  by  an  asymptotic  relaxation 
to  the  final  combustion  temperature.  As  a  representative  equation  for  exploring  the  implementation  of  the  POD- 
Galerkin  method  towards  the  development  of  a  reduced  order  model  for  combustion  applications,  we  consider  the 
one-dimensional  reaction-advection  scalar  equation. 


dt  dx 


=  K 


T{xj) 


(1) 


In  this  equation,  T  is  a  temperature-like  parameter  (hereafter  called  the  temperature)  and  X  represents  the  velocity  with 
which  the  temperature  is  convected  through  the  domain.  The  terms  on  the  right  hand  side  represent  a  competition 
between  a  source  term  that  is  dominant  at  lower  temperatures  and  a  sink  that  dominates  at  higher  temperatures.  The 
resulting  solution  consists  of  a  rapid  increase  in  temperature  from  an  initial  temperature,  =  T  (O) ,  followed  by  an 

asymptotic  approach  to  a  final  temperature,  =r(oo),  a  process  which  is  characteristic  of  a  typical  combustion 

problem.  The  parameter,  Tt,  in  the  second  term  controls  where  the  shift  between  source  and  sink  occurs  within  the 
flame  and  is  discussed  later. 


The  function,  t  (x)  ,  in  Eq.  1  is  the  solution  to  the  steady  version  of  Eq.  1  and  can  be  obtained  in  analytical  form 
as. 


f{x) 


'l  +  C,exp[fe] 


(2) 


The  constant. 


is  given  in  terms  of  the  initial  and  final  temperatures. 


To  avoid  extending  the  computational  domain  over  an  infinite  distance,  we  extend  the  domain  only  to  a  location 
where  the  temperature  has  very  nearly  reached  7}.  Specifically,  the  length  of  the  computational  domain  is  chosen  such 
that  at  X  =  L,  the  temperature  at  the  exit  of  the  domain  is  f{^L)  =  where  ^  «  1.  A  typical  value  used  for 

the  present  calculations  is  1%.  Having  determined  ^and  L,  the  scaling  constant,  ku  is  then  given  as: 


K  = 


sQ 


1  7 


(3) 


The  final  parameter,  Tt,  in  Eq.  (1)  is  a  constant  to  control  the  location  of  the  shift  between  source  and  sink  in  the 
linearized  version  of  the  equation  presented  below. 

Eig.  2  gives  a  representative  family  of  solutions  to  Eq.  (2)  for  an  initial  temperature  of  300  and  two  final 
temperatures,  2500  and  3000  for  three  different  values  of  s,  0.5%,  1%  and  2%,  for  each  final  temperature.  The 
resulting  steady  state  profiles  have  the  general  character  of  an  Arrhenius  term  with  a  rapid  increase  followed  by  an 
asymptotic  approach  to  the  final  temperature.  Values  of  ^  of  1%  imply  that  the  temperature  reaches  to  within  1%  of 
the  final  temperature  in  the  distance  x  =  L. 

The  following  sections  present  solutions  to  the  complete  non-linear  version  as  well  as  a  linearized  perturbation 
version  of  Eq.  (1).  The  perturbed  version  of  Eq.  (1)  is: 


dT'{x,t)  dXT'{x,t) 


dt 


dx 


i+ZL_2ZVlV'(x,f) 


(4) 
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where  T{xj)  =  r(jc)  +  r'(jc,r)  and  TIT  « 1 .  For  convenience,  we  represent  the  forcing  function  on  the  right  hand 
side  of  the  linearized  equation  as  the  function,  G: 


G  =  k, 


T{x)] 


(5) 


Note  that  G  changes  from  a  source  to  a  sink  at  the  particular  location  where  T  (v)  =  (^7^  +7^, )  /  2 .  When  Tt  is  zero 

this  transition  location  coincides  with  the  inflection  point  in  T  (x) .  The  transition  location  can  be  moved  to  the  right 
by  choosing  positive  values  of  Tj. 


Fig.  2  Characteristics  of  the  steady  state  solution,  t  (x)  versus  x,  for  two  final  temperatures  and  three  values 
of  s.  Solid  lines,  7}  =  3000K;  dashed  lines,  7}  =  2500K  {Tt=  0). 


B.  Model  Reduction  using  the  POD-Galerkin  Method 

A  flux-splitting,  finite  volume  scheme  with  uniform  cell  size  is  used  to  obtain  ‘high  fidelity’  solutions  of  both  the 
linear  and  nonlinear  versions  of  the  model  equation.  The  numerical  solutions  are  obtained  by  means  of  a  second- 
order  upwind  scheme  in  space  and  a  4*-order  Runge-Kutta  scheme  in  time.  For  a  typical  computation,  the  numerical 
solution  is  started  from  the  analytical  steady  state  solution  and  is  computed  time- accurately  for  a  sufficient  number 
of  steps  with  a  constant  flux,  specified  at  the  upstream  boundary  condition,  until  a  steady  numerical  solution  (which 
is  near,  but  not  exactly  equal  to,  the  analytical  solution)  is  reached.  The  upstream  boundary  condition  is  then 
switched  to  a  temporally  varying,  periodic  condition,  AT(0,t)  =  Q(t),  and  the  time-accurate  computation  is  continued 
until  stationary  conditions  have  been  reached.  The  solution  is  tabulated  and  stored  at  every  grid  point  for  a  large 
number  of  time  steps  (typically  on  the  order  of  10,000)  thereby  generating  a  rectangular  matrix  that  can  be  used  a 
data  base  for  fitting  an  eigen-basis  by  means  of  the  POD  procedure.  The  POD  basis  vectors  are  then  applied  to  the 
governing  linear  or  non-linear  PDF  to  derive  the  reduced-order  ODE  formulation.  The  procedure  for  deriving  the 
model  reductions  for  the  linear  and  non-linear  model  equations  is  described  in  the  next  two  sub- sections. 


a.  Model  reduction  of  linearized  equation 

The  linearized  version  of  Eq.  (1)  is  first  used  to  study  the  implementation  of  the  POD-Galerkin  method.  POD  is 
used  to  represent  variable  r'(x,^)  in  eigen-basis  format, 

5 
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(6), 


^  P  ^  P 

T'{x,t)  *  ^a„(f)(T„(4,(x)  =  ^a„(f)4(x) 

n=\  n=l 

where  cr„  is  the  singular  value  of  the  n*  POD  mode,  an(t)  is  the  n*  POD  temporal  mode  and  the  n*  eigen-basis, 
is  an  orthonormal  function. 


1,  if  k  =  n 
0,  otherwise 


(7). 


Singular  value  decomposition  (SVD)  is  used  to  obtain  the  eigen-basis  (pn(x).  Details  of  SVD  applications  for  obtaining 
the  POD  eigen-basis  can  be  found  in  Ref  25. 

Upon  obtaining  the  eigen-basis,  the  target  governing  equation  is  projected  onto  the  eigen-basis  (pk(x)  throughout 
the  whole  computational  domain,  X, 


j(4W 


dr{x,t)  _^dAr{x,t) 


dt 


dx 


X  V 

Substituting  the  POD  expansion,  Eq.  (6)  into  Eq.  (8)  gives. 


f  / 

=  j  (4  (•^) 

^I 

X 

V  V 

Tt 

1  +  ^-2 


Tf  7>  , 


T'{x,t) 


dx 


(8). 


X\  lA(x)^n(x)dx  +  \Aix) 

n=l  yx  j  n=l  yx 


dX(l)^  (x) 
dx 


dx 


n=l 


f 

r  f 

\A(x) 

K 

K  V 

1  +  ^-2 


f{x) 


J 


A  y 

\ 

dx 

)  J 

y 

(9). 


««(0 


Recognizing  that  the  eigen-basis  is  given  in  tabulated  form  requires  that  these  integrals  be  expressed  as  summations 
over  the  number  of  discrete  points  in  the  eigen-basis.  Eor  an  eigen-basis  tabulated  at  the  NI  cell  centers  of  the  finite 
volume  solution,  a  numerical  quadrature  on  the  remaining  integrals  gives. 


NI  A 

n=l  V  ^=1 


da„  (r) 


dt 


■X  Xaj 

n=l  V  ^=1 
(  NI 


dXtpn 

dx 


\ 

Ax^ 


^nit) 


f 

f  TV 

—  \ 
Ti 

\ 

h 

V 

-2  — 

Aj 

Axi 

Tf 

V  J 

^/J 

J 

J 

\  (0 


(10). 


Eor  notational  simplicity,  we  have  here  indicated  the  numerical  representation  of  the  flux  derivatives  symbolically. 
Both  central  and  upwind  differences  have  been  used  in  the  applications.  The  upstream  condition  is  obtained  from  the 

^P 

specified  boundary  condition,  (0^n,i/2  ^  ^(0  ’  thereby  providing  a  natural  way  for  incorporating  the  boundary 

n=l 

condition.  The  flux  on  the  downstream  boundary  (when  needed)  is  obtained  by  extrapolation,  taking  into  account  the 
hyperbolic  character  of  the  equation.  By  separating  the  contributions  of  boundary  condition  explicitly  as  a  function 
Fk(Q(t)),  we  can  get. 


y  NI 

n=l  V  1=1  J 

f 


dan  (t) 


dt 


'Yj  -Yaj 


n=l  \  i 


=1 


dx 


NI 
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J 

AXi 
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(11). 


an{t)  =  F,{a{t)) 
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Here  it  should  be  noted  that  operator  notation,  — ^  has  been  used  to  represent  the  gradient  term  left  (uncorrelated  to 

Sx 

boundary  condition)  after  separating  the  contributions  of  boundary  term  explicitly  in  Eq.  (10)  and  for  most  of  the 

^(  )  ^(  ) 

interior  cells  — ^  ^  but  can  be  different  at  cells  near  upstream/downstream  boundaries  and  dependent  of 

Sx  dx 

discretization  schemes  used. 

These  steps  reduce  the  PDE  in  Eqs.  (1)  or  (4)  to  the  system  of  ODE’s, 

d2i{t) 


M 


dt 


NI 


^  (O)  =  ^^4 =  O) Ax-  (initial  conditions) 

i=l 

where  =  a2{t)  •••  (^)]  . 

h(^)  =  (^)  [t)  •  •  •  with  contributions  from  boundary  conditions,  hk{t)  =  FdSl{t)), 

M  is  identity  matrix,  (M  =  I),  due  to  the  orthogonality  of  POD  eigen-basis. 

The  elements  in  matrix  L  are  determined  by  the  numerical  quadrature  related  to  an(t)  in  Eq.  (11), 

NI  /  OT  j  \  NI 


(12) 


^kn  ^^j^k,i 


i=l 


Sx 


i=l 


^  ^  T  f^  ^ 

V  )  J 


K 

V  V 


AX;  . 


b.  Model  reduction  of  nonlinear  equation 

The  application  of  the  POD-Galerkin  method  to  the  nonlinear  model  equation  (Eq.  (1))  is  very  similar.  The  POD 
expansion  is  used  to  approximate  the  fluctuating  quantity  by  extracting  the  time-averaged  value  from  the  numerical 
solution, 

T{x,t)  =  T^^^{x)  +  T'{x,t)«T^^^  (q  + =  r^ve(x)  + (13), 

n=l  n=\ 


where  rave(x)  is  the  time-averaged  value  of  T(x,t).  Note  that  T^yeix)  is  the  average  of  the  non-linear  solution  and,  in 
general,  is  not  equal  to  the  analytical  function,  T  (x) . 


The  same  procedure  can  be  followed  from  Eq.  (8)  to  Eq.  (1 1)  to  apply  model  reduction, 

dan  (0 


Np  r  NI 


n=l  t 1=1 


dt 


Np  \  NI 
n=l  i=l 


i  5x 


h 


m=l  n=\  1  1=1 


y  '^f  j 


1=1 

l^n  (^)^m  (^) 


^  f  T  f.'^ 

1  +  ^-2-^ 

V  V  y  j 


(14). 


1=1 


dx 


rk,i 


f  rj.  rj^2  ^ 

kT  -d—(T  —T\—h 

N^ave,i^N^  (^ave,i  ^i)  N  rj. 


)  1 

y. 

)  \ 

The  reduced  lower-order  ODE  is. 
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(15) 


M  - L^a(0  +  n(^a(f))  =  h(f)-b 

NI 

a,{0)  =  Y,<t>ki’^i)T'Wt  =  O)  Ax-  (initial  conditions ) 


where  a(^)  =  |^ai(^)  a^{t)  •••  (^)]  . 

h(^)  =  {t)  {t)  •  •  •  with  contributions  from  boundary  conditions,  hk{t)  =  Fk{0.{t)), 

b  =  Z?2  •  •  •  j  with  contributions  from  time-average  value  TaveC-x), 


h  ^l7;ve,,+*i:^(7;ve,,-7])-*l^  Ax, 

If  If 


Ln  is  similar  to  matrix  L  in  Eq.  (12)  but  the  time  averaged  value,  Ts^y^ix)  should  be  used  to  replace  t  (x)  , 

n(a(0)  is  a  quadratic  term  coming  from  the  nonlinear  term  of  the  model  equation, 

(  m  ( k6  .6  ^ 

=  ‘  "j  Ax;L„(f)a„(f) 

m=l  n=\  1  i=\  \  f  J  \ 


III.  Overview  of  Test  Problem 

In  selecting  the  parameters  in  the  model  equation,  and  the  computational  domain,  we  have  chosen  values  that  are 
representative  of  a  typical  combustion  problem.  The  computational  domain  goes  from  0  to  0.5  m  L  =  0.5m) 
(representative  of  a  practical  combustor  length)  with  1,000  uniform  cells.  The  values  of  model  constants  (To,  Tf,  ^and 
Tt)  in  the  right-hand-sides  of  Eq.  (1)  and  Eq.  (4)  have  been  selected  as  To  =300K  and  7}  =3000K  corresponding  to 
cold  incoming  reactants  and  a  nominal  flame  temperature.  The  parameter,  s,  has  been  set  to  0.01,  allowing  the 
temperature  to  reach  within  1  %  of  the  final  temperature  inside  the  computational  domain.  These  three  quantities  yield 
the  mean  profile  shown  in  Eig.  3.  Einally,  the  transition  temperature,  Tt,  has  been  set  to  1500K,  and  the  corresponding 
variation  of  the  function,  G,  for  the  linear  version  (Eq.  (5))  is  also  shown  in  Eig.  4.  This  value  of  Tt  has  been  selected 
so  that  the  location  where  G  transitions  from  driving,  G  >  0,  to  damping,  G  <  0,  in  the  linear  equation  occurs  at  x  = 
0.24m  (approximately  the  middle  of  the  computational  domain)  as  indicated  by  the  dashed  line.  The  same  value  of  Tr 
is  also  used  for  the  non-linear  equation  resulting  in  a  qualitatively  similar  behavior.  The  convection  velocity,  X  has 
been  set  as  50  m/s,  a  nominal  speed  in  a  practical  combustor.  The  time  step  used  for  numerical  solution  is  10'^  s  and 
the  solutions  are  output  every  10  time  steps  (10'^  s). 

Eluctuating  conditions  inside  the  computational  domain  were  obtained  by  specifying  a  periodic  upstream  boundary 
condition,  Q(0,  of  the  form. 


Jsin(2;r/,?) 

k=\ _ 

max  ^sin(2;r/^/) 


XT(^Xy^,t)- XTq  1-h^^ 


where represents  the  frequency  of  the  N  modes  that  are  used  in  the  boundary  condition.  The  number  of  modes,  N, 
used  in  forcing  ranges  from  a  single  sinusoidal  mode,  N  =  1,  to  a  total  of  eight  modes,  /V  =  8.  The  peak  value  of  the 


multi-frequency  sinusoidal  signal,  max  j^^sin  (2;r/^^)J  is  used  to  normalize  the  fluctuating  coefficient  to  make  sure 
the  total  fluctuating  component  is  Sf  of  the  reference  value  XTo. 
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Mean  Temperature 
Forcing  Coefficient  G 


Fig.  3  Mean  temperature  profile  and  forcing  coefficient  G  of  test  case  (G  >  0,  driving;  G  <  0,  damping). 


In  the  results  reported  in  the  following  section,  we  first  start  with  solutions  from  the  linear  equation,  Eq.  (4)  and 
then  consider  the  non-linear  Eq.  (1).  The  initial  examples  from  the  linear  equation  provide  insight  without  the 
complexities  of  non-linearities  while  the  non-linear  examples  allow  us  to  concentrate  on  any  new  behavior  related  to 
the  non-linearity.  Eor  both  linear  and  non-linear  examples,  the  model  reduction  studies  follow  three  distinct  steps. 
Eirst,  we  verify  the  POD  expansion  with  the  original  CED  solution  to  ensure  that  it  accurately  reproduces  the  numerical 
solution.  Second,  we  compare  the  ROM  solution  against  the  original  CED  solution  that  was  used  to  derive  the  POD 
basis  vectors.  This  allows  us  to  ascertain  how  well  the  ROM  can  reproduce  the  CED  solution.  Third,  we  extend  the 
ROM  to  conditions  different  from  the  original  training  CED  solution  and  then  use  new  CED  calculations  to  verify  the 
extended  ROM  solutions.  The  final  test  is  important  because  this  showcases  how  the  ROM  will  actually  be  used  within 
a  design  study,  i.e.,  to  investigate  the  effects  of  changes  in  geometric  configurations  and/or  operating  conditions. 


IV.  Model  Reduction  Studies  using  POD  Eigen-basis 

Studies  of  model  reduction  techniques  are  covered  in  this  section  in  two  parts: 

1 .  Validation  of  ROM  using  single  frequency  forcing  (linearized  model  equation) 

2.  Validation  of  ROM  using  broadband  frequency  forcing  (both  linearized  and  nonlinear  model  equation) 

A.  ROM  validation  using  single  frequency  forcing 

As  a  first  step,  the  linear  model  equation  (Eq.  (4))  is  solved  by  specifying  a  single  frequency  forcing  at  the  upstream 
boundary.  The  resulting  numerical  solutions  are  then  used  to  calculate  the  POD  eigen-basis  upon  which  the  ROM  is 
based.  The  boundary  condition  is  perturbed  at  a  single  frequency  500Hz  with  Sf=  0.01  (Eq.(16))  to  obtain  numerical 
solutions  by  solving  the  linearized  model  equation.  The  results  at  four  locations  are  shown  in  Pig.  4.  It  is  evident  that 
the  fluctuating  amplitude  increases  as  it  goes  downstream  till  the  inflection  point  shown  in  Pig.  3  (i.e.,  the  source  term 
causes  growth)  and  once  it  passes  that  point,  the  amplitude  decreases  (i.e.,  the  source  terms  causes  damping). 

In  the  following  sub-section,  we  first  verify  the  solutions  obtained  from  the  POD  eigen-basis.  Pollowing  this,  the 
original  numerical  solutions  are  then  used  to  validate  the  reconstructed  solutions  from  the  ROM.  In  the  next  sub¬ 
section,  we  evaluate  the  robustness  of  the  ROM  procedure. 

a.  Verification  of  ROM  procedure 
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X  =  0.25  m  x  =  0.375  m 

Fig.  4  CFD  results  from  linearized  model  equation  with  500Hz  boundary  perturbation  at  four  different 

locations. 


To  investigate  the  quality  of  the  POD  eigen-basis,  the  power  in  any  individual  mode,  n,  (herein  called  the  Mode 
Power)  and  the  cumulative  sum  of  the  power  of  all  modes  up  to  Np  (the  Cumulative  Power)  are  defined  based  on  the 
6-^  value  in  Eq.  (6)  and  given  as  below. 


Mode  Power 


^p, total 

z 


n=\ 


xl00% 


Cumulative  Power  ( - x  100% 

V  P  J  Nppotal 


(17). 


where  Np  is  the  number  of  POD  modes  used  in  reconstructing  the  solutions.  These  power  definitions  are  commonly 
used  in  POD  analyses.  From  Fig.  5(a),  both  the  individual  POD  mode  power  and  the  cumulative  power  indicate  that 
the  first  two  POD  modes  include  nearly  100%  of  the  information  in  the  numerical  solution,  while  the  energy  content 
in  all  the  higher  modes  is  ten  orders  of  magnitude  smaller. 

In  addition  to  denoting  the  power  in  the  various  POD  modes,  we  also  introduce  an  L2-norm  based  upon  the  error 
between  the  eigen-basis  expansion  and  the  original  CFD  solution.  This  L2  error  provides  a  more  sensitive  assessment 
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of  how  well  the  reconstructed  solutions  represent  the  original  CFD  result.  The  L2  error  at  any  location  Xi  is  defined 


as, 


L2-norn[i 


K)=- 


jT'\x.,t)dt 


(18) 


As  shown  in  Fig.  5b,  including  only  the  first  two  POD  modes  results  in  an  L2  error  of  10'^.  Interestingly,  adding 
more  modes  does  not  decrease  the  L2  error,  indicating  that  the  energy  content  in  the  higher  modes  for  this  single 
wave  forcing  problem  is  indeed  very  small.  As  shown  later,  multiple -mode  forcing  requires  substantially  more 
modes  to  reach  this  level  of  error. 


It  is  important  to  note  that  in  addition  to  using  the  L2-norm  defined  in  Eq.  (18)  to  compare  the  eigen-basis 
expansion  with  the  original  solution,  it  can  also  be  used  to  compare  the  ROM  solution  with  the  CFD  solution.  To 
assess  the  accuracy  of  the  ROM  solution,  the  temporal  modes.  Unit),  in  Eq.  (18)  are  taken  from  the  ROM  solution 
(Eq.  (12)),  instead  of  from  the  POD  expansion  (Eq.  (6)).  In  either  case,  the  spatial  modes,  (l)n{xi)  are  taken  from  the 
POD  eigen-basis  (Eq.  (6)). 


x  =  0m  x  =  0.125m 


Mode  #  Mode  # 


(a)  (b) 

Fig.  5  Validations  of  POD  eigen-basis  from  linearized  model  equation  with  500Hz  boundary  perturbation  (a) 
POD  mode  power  &  cumulative  power  versus  mode  number  (b)  L2-norm  of  the  reconstructed  solution  from 

POD  at  four  locations. 

To  characterize  the  reduced  order  model,  the  eigenvalues  of  stiffness  matrix  L  (cr  +  72;^  where  cr  is  the  real  part 
indicating  growth  rate,  7  is  the  square  root  of  minus  one  and  2;^ is  the  imaginary  part  indicating  the  frequency,/)  in 
Eq.  (12)  are  calculated  and  shown  in  Fig.  6  using  a  2“^-order  upwind  scheme  to  approximate  eigen-basis  gradient, 
which  is  consistent  with  the  CFD  solutions.  With  the  increase  of  POD  modes  in  building  the  ROM,  the  responses  at 
the  target  forcing  frequency  (500Hz,  Fig.  6(b))  are  consistent  (cr=  -1.9)  while  as  more  POD  modes  are  included  in  the 
ROM,  responses  at  other  frequencies  show  up  but  they  are  situated  far  from  the  zero-axis  (cr=  0)  and  more  importantly 
sit  on  the  left  half  of  the  plane.  This  indicates  they  are  all  highly-damped/stable  modes  that  have  little  influence  on  the 
ROM  solution. 

ROM  calculations  based  upon  using  central  differences  to  approximate  the  eigen-basis  gradient  in  Eq.  (10)  are 
similar,  but  found  to  be  inferior  to  upwind  differencing.  Similar  responses  at  the  forcing  frequency  are  found  but 
unlike  in  the  case  using  2“^-order  upwind,  as  more  POD  modes  are  included,  neutrally  stable  (near  the  zero  axis)  or 
even  unstable  modes  (right  half  of  the  plane)  at  other  frequencies  show  up.  It  is  believed  that  this  difficulty  is  not  a 
characteristic  of  central  differencing,  but  rather  arises  from  using  inconsistent  differencing  between  the  CFD  and  ROM 
simulations.  The  effect  of  these  problematic  modes  on  the  ROM  solution  is  covered  at  the  end  of  this  section. 
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Fig.  6  Eigenvalues  of  the  stiffness  matrix  L  in  ROM  (500Hz  forcing)  (a)  overview  of  all  eigenvalues  and  (b) 
zoomed-in  region  with  dominant  eigenvalues  (cr>  0  means  unstable  modes  and  cr<  0  means  stable  modes). 


The  reconstructed  solution  using  ROM  with  one,  two  and  five  modes  is  compared  with  the  CFD  solution  near  the 
inflection  point  of  the  domain  in  Fig.  7.  Here  the  ROM  is  built  using  2”^-order  upwind  to  approximate  the  eigen-basis 
gradient.  The  solutions  are  reconstructed  in  the  same  way  defined  in  Eq.  (6)  but  the  POD  temporal  modes,  an(t)  are 
calculated  from  the  ROM  (Eq.  (12))  with  the  POD  eigen-basis,  (pn(x)  calculated  from  the  CFD  solutions.  Consistent 
with  the  findings  in  the  POD  validation  (Fig.  5),  after  including  more  than  2  POD  modes,  the  ROM  solution  shows 
very  good  agreement  with  CFD  solution. 


Fig.  7  Comparisons  between  the  reconstructed  ROM  solutions  and  the  CFD  solutions  at  jc  =  0.25  m  with  At  = 

le-5  s. 


To  provide  a  more  quantitative  assessment  of  how  the  ROM  predictions  compare  with  the  CFD  results  we  present 
in  Fig.  8  the  L2-norms  of  the  ROM  solution  at  the  mid-point  of  the  computational  domain.  For  reference,  the  L2  norm 
for  the  POD  expansion,  originally  given  in  Fig.  5,  is  also  included.  Fig.  8(a)  shows  the  L2  norms  for  ROM  time  steps 
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of  At  =  le-5  and  le-6  s.  (Recall  that  the  original  numerical  solutions  were  simulated  with  At  =  le-6  s  while  the 
solutions  are  output  every  10  time  steps  so  the  POD  eigen-basis  are  built  on  numerical  solutions  with  At  =  le-5  s.)  For 
both  values  of  the  ROM  time  step,  the  convergence  of  the  ROM  L2-norm  is  consistent  with  that  from  the  POD  norm. 
The  error  drops  dramatically  from  the  first  to  the  second  mode,  but  there  is  no  further  improvement  when  additional 
modes  are  added.  There  is,  however,  a  difference  in  level.  For  a  time  step  of  le-5  s,  the  L2  error  for  two  or  more 
modes  is  nominally  10'^  while  for  At  =  10'^  the  L2  error  is  10'^.  Putting  this  in  perspective,  the  L2  error  for  the  POD 
expansion  was  10'^.  These  results  suggest  the  ROM  solution  is  second-order  in  time  (as  would  be  expected  from  the 
time  discretization). 

The  L2  norm  of  the  error  in  the  ROM  solution  also  depends  upon  the  method  used  for  discretizing  the  function 
derivatives  in  Eq.  (10)  as  can  be  seen  by  comparing  Fig.  8(b)  with  Fig.  8(a).  The  results  in  Fig.  8(a)  are  based  upon 
using  a  second-order  accurate  discretization  which  is  consistent  with  that  used  in  the  CFD  solution,  whereas  those  in 
Fig.  8b  are  for  central  difference.  The  error  convergence  in  the  central  difference  case  is  qualitatively  for  both  time- 
step  sizes,  but  the  levels  of  convergence  reached  are  very  different.  The  error  in  the  central  difference  calculation  stalls 
at  10'^  for  time  steps  of  both  10'^  and  10'^  s.  The  inconsistency  between  ROM  solution  using  central  difference  and 
CFD  solution  can  be  attributed  to  the  problematic  neutrally  stable/unstable  eigen-modes  mentioned  above.  Therefore, 
to  build  up  a  ROM  that  is  able  to  give  reasonable  and  comparable  predictions  of  the  physics  of  interest,  it  appears  that 
the  numerical  scheme  used  in  building  the  ROM  should  be  consistent  with  that  used  in  calculating  the  CFD  solutions. 


(a)  2”^-order  upwind  (b)  central  difference 

Fig.  8  L2-norm  of  the  reconstructed  solutions  from  ROM  at  x  =  0.25  m  using  different  numerical  scheme  in 
reconstructing  eigen-basis  at  cell  faces  with  (a)  2”*^-order  upwind  (b)  central  difference. 


b.  Effects  of  time  step  in  the  ROM 

To  further  investigate  the  influence  of  the  time  step  size  on  the  ROM  solutions,  the  ROM  is  built  from  CFD 
solutions  under  a  much  lower  frequency  of  excitation,  i.e.,  20Hz  with  Sf=  0.01.  The  second-order  upwind  scheme  is 
used  in  ROM  to  approximate  the  eigen-basis  gradient.  Different  time  steps  are  used  to  calculate  ROM  solutions.  Low- 
frequency  excitation  (20Hz)  is  used  to  make  sure  that  the  time  steps  used  can  provide  sufficient  sample  rate  for  the 
target  frequency.  For  a  20Hz  wave,  a  minimum  sample  rate  of  40Hz  or  0.025  s  is  needed  to  avoid  aliasing.  Similar 
POD  validation  procedures  have  been  followed  as  in  Fig.  5  and  it  has  again  been  found  that  two  POD  modes  are 
sufficient  to  capture  the  CFD  solution. 

In  this  study,  two  sets  of  POD  eigen-basis  have  been  calculated  using  two  different  sample  rates  (sampled  every 

10  time  steps,  Atpod  =  le-5  s  and  100  time  steps,  Atpod  =  le-4  s)  of  CFD  solutions  respectively  to  build  the  ROM’s.  The 
eigenvalues  of  the  ROM  stiffness  matrix  L  are  shown  in  Fig.  9  and  Fig.  10.  The  predictions  of  the  responses  at  target 
frequency  (20Hz)  are  consistent  and  similar  to  the  previous  case  (Fig.  6).  Responses  at  frequencies  other  than  the 
driving  frequency  are  highly  damped.  Though  the  two  sets  of  ROM  are  built  from  different  POD  eigen-basis,  the 
information  delivered  from  the  ROM’s  is  similar. 

The  comparisons  between  the  ROM  solutions  with  different  time  steps  and  the  CFD  solutions  are  shown  in  Fig. 

1 1  with  Atpod  =  le-5  s.  CFD  solutions  are  compared  with  the  ROM  solution  using  two  and  four  POD  eigen -basis  at  x 
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=  0.25  m.  In  Fig.  11(a)  with  relatively  higher  temporal  resolution  (At  =  le-6,  le-5,  le-4  and  le-3  s),  the  predictions 
from  the  ROM  by  only  including  two  POD  eigen-basis  follow  the  numerical  solution  well  qualitatively.  As  the 
temporal  resolution  gets  poorer,  the  ROM  is  still  able  to  capture  the  same  phase  as  the  numerical  solution  but 
deviations  in  amplitude  predictions  show  up  (At  =  5e-3  and  le-3  s).  As  the  temporal  resolution  gets  close  to  the 
required  sample  rate  (At  =  2e-2  s),  the  prediction  from  the  ROM  fails  to  follow  the  CFD  solution.  By  using  a  ROM 
built  from  4  POD  eigen-basis  (Fig.  1 1(b)),  the  ROM  solutions  with  higher  temporal  resolution  ((At  =  le-6  and  le-5  s) 
reach  good  agreement  with  CFD  solution  while  by  using  larger  time  steps,  the  ROM  solutions  fail  to  predict  the  correct 
behavior. 
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Fig.  9  Eigenvalues  of  the  stiffness  matrix  L  in  ROM  (20Hz  forcing)  built  from  POD  eigen-basis  calculated 
using  CFD  results  sampled  every  10  time  steps  (Atpod  =  le-5  s)  (a)  overview  of  all  eigenvalues  and  (b)  zoomed- 
in  region  with  dominant  eigenvalues  (  o'>  0  means  unstable  modes  and  C7<  0  means  stable  modes). 
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Fig.  10  Eigenvalues  of  the  stiffness  matrix  L  in  ROM  (20Hz  forcing)  built  from  POD  eigen-basis  calculated 
using  CFD  results  sampled  every  100  time  steps  (Atpod  =  le-4  s)  (a)  overview  of  all  eigenvalues  and  (b) 
zoomed-in  region  with  dominant  eigenvalues  ((T>  0  means  unstable  modes  and  cr<  0  means  stable  modes). 
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(a)  Using  2  POD  eigen-basis 


(b)  Using  4  POD  eigen-basis 


Fig.  11  Comparisons  between  the  ROM  solutions  and  the  CFD  solution  at  jc  =  0.25  m  with  Atpod  =  le-5  s  using 
(a)  2  POD  eigen-basis  (b)  4  POD  eigen-basis,  (c)  L2-norm  of  the  ROM  solutions  at  jc  =  0.25  m  with  different 

time  steps. 


This  can  be  further  confirmed  quantitatively  by  looking  at  the  L2-norm  of  ROM  solutions  at  x  =  0.25  m  with 
different  time  steps  in  Fig.  12.  For  time-steps  larger  than  le-4  s,  the  ROM’s  built  with  more  than  two  POD  eigen- 
bases  diverge  and  cannot  be  properly  used  with  the  POD  eigen-basis  calculated  from  higher  sample  rate  (Fig.  12(a)). 
Even  when  a  lower  sample  rate  (Atpod  =  le-4  s)  is  used  to  calculate  the  POD  eigen-basis  from  the  CFD  solutions,  the 
ROM  solutions  still  diverge  in  the  same  fashion  (Fig.  12(b)),  which  might  indicate  that  the  temporal  resolution  of  the 
POD  eigen -basis  construction  is  less  critical  in  determining  the  ROM  predictions.  However,  the  failure  of  the  ROM 
built  with  more  POD  eigen-bases  can  be  attributed  to  a  lack  of  temporal  resolution  to  describe  those  highly  damped 
modes  as  shown  in  Fig.  9  and  Fig.  10.  With  both  sample  rates,  the  growth  rates,  a  of  those  modes  are  predicted  to  be 
less  than  -1,00,000,  which  indicates  a  decay  rate  of  1/1,00,000  =  le-5  s.  By  including  more  than  2  POD  eigen-bases 
to  build  the  ROM,  higher  temporal  resolution  is  required  to  capture  those  large  decay  rates . 
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(a)  Atpod  =  le-5  s  (b)  Atpod  =  le-4  s 

Fig.  12  Comparisons  between  the  ROM  solutions  and  the  CFD  solution  at  x  =  0.25  m  with  (a)  Atpod  = 
le-5  s  (b)  Atpod  =  le-4  s  at  x  =  0.25  m  with  different  time  steps. 
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In  summary  of  this  section,  reasonable  predictions  are  obtained  even  when  larger  time  steps  are  used  for  the  ROM 
when  it  is  built  using  two  POD  eigen-bases.  However,  when  more  than  two  POD  eigen-bases  are  used  then  higher 
temporal  resolution  is  required,  ostensibly  to  capture  the  large  decay  rate  at  other  frequencies. 


B.  ROM  Implementations  using  broadband  forcing 

The  studies  above  are  all  focused  on  the  validation  of  the  ROM.  In  this  section,  the  implementation  of  the  ROM 
for  predicting  the  response  at  other  forcing  frequencies  is  investigated.  To  do  that,  a  numerical  solution  is  obtained 
with  broadband  forcing  (200-l,000Hz  with  lOOHz  frequency  increments,  Af,  and  Sf=  0.01)  at  the  boundary  so  that  the 
ROM  built  is  expected  to  be  capable  of  capturing  the  response  at  those  frequencies.  A  set  of  separate  CFD  solutions 
has  been  calculated  with  different  forcing  under  a  range  of  frequencies  as  reference  cases  to  test  the  capabilities  of  the 
derived  ROM  (summarized  in  Table  1).  Both  linear  and  non-linear  cases  are  considered. 


Case  Information 

Forcing  frequencies,  Hz 

Case  used  to  build  ROM 

Broadband  (200-l,000Hz  with  Af  =  lOOHz) 

Test  cases  for  ROM  to  predict 

100 

200 

500 

1000 

200  &  500 

200,  500  &  1000 

200,  500,  800  &  1000 

Table  1.  Summary  of  cases  used  in  studies  of  ROM  implementations. 


a.  ROM  trained  from  linear  model  equation 

The  representative  CFD  results  for  broadband  forcing  are  shown  in  Fig.  13  at  four  locations.  The  multiple  modes 
result  in  a  periodic  fluctuation  that  is  much  different  from  the  single  frequency  forcing  in  Fig.  4,  but  the  overall 
response  is  similar:  the  fluctuating  amplitude  grows  as  it  goes  downstream  till  the  inflection  point  and  decays  after  it 
passes  that  point. 

The  POD  mode  power  and  the  cumulative  power  for  this  multi-frequency-forced  case  are  shown  in  Fig.  14(a).  It 
is  particularly  interesting  to  compare  these  results  with  those  from  the  single  frequency  case  in  Fig.  5.  In  both  the 
present  multi-frequency  case  and  the  single  mode  case,  there  is  a  grouping  of  modes  with  a  substantial  amount  of 
energy  followed  by  a  sudden  decrease  to  negligible  energy  in  the  remaining  higher  modes.  For  the  single  frequency 
case,  the  first  two  modes  contained  nearly  100%  of  the  energy  with  about  50%  in  each  mode.  For  this  present  case, 
the  first  18  modes  contain  nearly  100%  of  the  energy,  with  some  10%  in  the  first  mode  followed  by  a  monotonic  decay 
until  the  18*  mode  contains  about  1%.  Modes  19  and  higher  all  contain  negligible  amounts  of  energy.  Recalling  that 
nine  modes  were  used  to  drive  the  present  case  while  a  single  mode  in  the  results  of  Fig.  5,  it  appears  that  two  POD 
modes  are  required  to  represent  each  driving  frequency.  The  cumulative  power  in  the  two  cases  is  also  consistent  with 
this  picture. 

This  conclusion  is  confirmed  by  the  L2-norm  of  the  reconstructed  POD  solutions  in  Fig.  14(b).  The  L2  error  remains 
large  until  a  total  of  1 8  modes  have  been  included,  at  which  time  the  L2  error  drops  dramatically.  Incorporating  more 
than  1 8  modes  provides  no  further  decrease  in  the  error  again  indicating  that  all  the  energy  is  represented  in  these  first 
18  modes. 
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Fig.  13  CFD  results  from  linearized  model  equation  with  broad  band  boundary  perturbation  at  four  different 

locations. 


POD  Mode  # 


X  =  Om 


x  =  0.125m 


(a)  (b) 

Fig.  14  (a)  POD  mode  power  &  cumulative  power  versus  mode  number  (b)  Li-norm  of  the  reconstructed 
solutions  using  POD  from  linearized  model  equation  with  broadband  forcing. 
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The  characteristics  and  validation  of  the  ROM  built  from  the  broadband  forcing  are  shown  in  Fig.  15.  The 
eigenvalues  of  the  stiffness  matrix  L  (Fig.  15(a)  and  (b))  indicate  that,  by  including  more  than  18  POD  modes,  the 
ROM  predictions  of  responses  at  the  target  forcing  frequencies  become  consistent.  Similar  to  the  observations  made 
in  the  single-frequency  cases,  there  are  some  highly-damped/stable  frequencies  sitting  in  the  left  plane.  The  L2-norm 
of  the  ROM  solution  at  v  =  0.25  m  (Fig.  15(c))  also  leads  to  similar  conclusions,  i.e.,  the  inclusion  of  eighteen  POD 
modes  yield  good  convergence  of  the  error.  Again,  for  the  smaller  time  step  (At  =  le-6  s),  the  accuracy  level  of  the 
ROM  solution  is  improved  and  closer  to  the  POD  solution. 

Finally,  as  a  test  to  see  if  the  ROM  trained  by  the  broadband  forcing  can  predict  the  correct  response  for  different 
forcing  frequencies,  we  try  the  cases  listed  in  Table  1.  The  different  frequencies  are  fed  to  the  ROM  by  modifying  the 
function  Q(0  in  Eq.  (12)  based  on  the  definition  in  Eq.(16).  The  L2-norm  of  the  ROM  solutions  are  shown  in  Eig.  16 
with  At  =  le-5  s,  indicating  reasonable  predictions  for  all  cases  except  the  100  Hz  forcing  case.  We  note  that  the  lOOHz 
forcing  is  outside  the  range  of  broadband  forcing  that  was  used  to  train  the  ROM  and  it  is  therefore  expected  that  this 
case  would  not  be  represented  correctly. 


(a) 


(b) 


Fig.  15  Characteristics  of  ROM  (broadband  forcing)  with  eigenvalues  of  the  stiffness  matrix  L  in  ROM  (a) 
overview  of  all  eigenvalues  and  (b)  zoomed-in  region  with  dominant  eigenvalues  (a  >  0  means  unstable  modes 
and  a  <  0  means  stable  modes)  (c)  Li-norm  of  the  ROM  solutions  at  jc  =  0.25  m . 
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Fig.  16  L2-norm  of  the  linear  ROM  solutions  at  four  locations  under  different  forcing  frequencies  with  At  = 

le-5  s. 


b.  ROM  trained  from  nonlinear  model  equation 

Similar  to  the  detailed  studies  using  the  linearized  governing  equation,  explorations  of  the  ROM  trained  from 
nonlinear  model  equation  are  shown  next.  As  before,  broadband  forcing  with  Sf  =0.01  is  used  to  generate  numerical 
solutions  for  construction  of  the  ROM.  Representative  CFD  results  at  different  locations  are  shown  in  Fig.  17.  For  the 
low  amplitude  forcing  that  is  used  here,  the  fluctuation  amplitudes  at  each  location  is  very  close  to  those  in  the  linear 
case  (Fig.  13).  Validation  of  the  POD  and  the  ROM  are  shown  in  Fig.  18.  A  big  drop  in  mode  power  can  be  observed 
for  POD  modes  higher  than  eighteen  and  the  cumulative  mode  power  reaches  close  to  100%  (Fig.  18(a)).  This  can  be 
interpreted  as  the  inclusion  of  complete  information  to  describe  the  linear  responses,  which  is  consistent  with  the 
finding  in  the  linear  case  (Fig.  14(a)).  However,  unlike  the  linear  case,  the  mode  power  is  observed  to  keep  decreasing 
as  we  go  to  even  higher  modes  and  saturates  around  70. 

The  validation  of  the  POD  and  ROM  reconstructions  shown  in  Fig.  18(b)  can  provide  more  quantitative 
understanding.  The  convergence  of  the  L2-norms  of  POD  and  ROM  is  consistent  with  the  POD  mode  &  cumulative 
power.  With  At  =  le-5  s,  the  L2-norms  of  ROM  saturate  after  40  POD  modes  and  converge  by  about  four  orders  of 
magnitude  while  by  using  At  =  le-6  s,  the  convergence  of  the  ROM  closely  follows  that  of  the  POD  solution  itself 
goes  down  by  nearly  eight  orders  of  magnitude. 
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Fig.  17  CFD  results  from  nonlinear  model  equation  with  broad  band  boundary  perturbation  at  four 

locations. 
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Fig.  18  (a)  POD  mode  power  &  cumulative  power  versus  mode  number  (b)  Li-norm  of  the  reconstructed 
solutions  using  POD  and  ROM  from  nonlinear  model  equation  with  broadband  forcing. 
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As  in  the  linear  case,  we  next  test  the  nonlinear  ROM  by  forcing  it  with  different  frequencies  as  shown  in  Table  1 
with  the  exception  that  we  leave  out  the  100  Hz  forcing.  The  L2-norm  of  nonlinear  ROM  solutions  are  shown  in  Fig. 
19  at  X  =  0.25  m.  It  is  noted  that  it  takes  more  POD  modes  to  reach  convergence  of  the  L2-norm  in  the  nonlinear  ROM 
than  in  the  linear  case.  In  Fig.  16,  only  18  modes  are  needed  to  attain  convergence  saturation  while  in  Fig.  19,  more 
than  40  modes  are  needed  to  reach  the  saturation  level.  This  might  be  attributed  to  the  nonlinearity  of  the  model 
equation.  In  addition,  unlike  the  linear  case  (Fig.  16),  where  predictions  of  responses  at  other  frequencies  can  reach 
better  accuracy  level  than  the  reference  case  (broadband  forcing),  in  the  nonlinear  case,  the  reference  case  achieves 
faster  convergence  and  a  better  accuracy  level.  Nevertheless,  we  conclude  that  the  nonlinear  ROM  is  capable  of 
predicting  responses  at  different  frequencies  to  a  reasonable  degree  of  accuracy. 
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Fig.  19  L2-norm  of  the  ROM  solution  at  jc  =  0.25  m  under  different  forcing  frequencies  with  At  =  le-5  s. 

For  the  nonlinear  cases,  we  note  that  the  forcing  amplitude  can  also  be  a  critical  variable  for  establishing  the 
capabilities  of  the  ROM.  To  test  that,  the  forcing  amplitudes  have  been  varied  with  Sf  =  0.01,  0.02,  0.04,  0.08  and 
0.16  for  the  same  broadband  frequencies  (sf  =  0.01  is  the  reference  case  to  build  ROM). 

The  middle  point  in  computation  domain  (x  =  0.25m)  has  been  picked  to  compare  the  ROM  solution  and  CFD 
solution  as  shown  in  Fig.  20.  The  Sf  =0.16  case  fails  to  follow  the  CFD  solution  while  the  other  three  cases  (sf  = 
0.02,  0.04  and  0.08)  capture  the  correct  phase  information  with  regards  to  the  CFD  solutions.  Sf  =  0.02  and  Sf  =  0.04 
reach  reasonable  agreement  in  the  fluctuating  amplitudes  also.  But,  ^  =  0.08  appears  to  over-predict  the  fluctuating 
amplitudes.  These  results  suggest  that  the  ROM  built  using  lower  amplitude  forcing  is  unable  to  properly  resolve  all 
of  the  non-linear  phenomena  encountered  at  higher  amplitudes. 

This  information  can  be  quantitatively  represented  in  the  L2-norms  of  the  nonlinear  ROM  solutions  for  different 
forcing  amplitudes  as  shown  in  Fig.  21.  Reasonable  accuracy  level  can  be  reached  by  including  more  than  40  POD 
modes  for  Sf  =  0.02,  0.04  and  more  than  50  modes  for  Sf  =  0.08  while  for  Sf  =0.16,  the  nonlinear  ROM  predictions 
are  observed  to  diverge  immediately  after  including  more  than  16  POD  modes.  One  interesting  finding  is  that  the 
cases  with  Sf  =  0.04  and  0.08  diverge  after  including  more  than  18  POD  modes  and  come  back  to  reasonable 
convergence  when  higher  POD  modes  are  included. 
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Fig.  20  Comparisons  between  the  ROM  solution  and  the  CFD  solution  at  x  =  0.25  m  using  different  forcing 

amplitude  with  At  =  le-5  s. 
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Fig.  21  L2-norm  of  the  reconstructed  solutions  from  nonlinear  ROM  at  jc  =  0.25  m  with  different  forcing 

amplitudes  with  At  =  le-5  s. 

Overall,  the  ROMs  built  from  the  linearized  and  nonlinear  model  equation  provide  useful  information  and,  more 
importantly,  are  capable  of  predictions  of  the  responses  for  different  forcing  amplitude  and  frequencies  without  having 
to  go  through  the  solution  of  the  full  PDE  model.  Thus,  the  ROM  procedure  can  provide  substantial  savings  in 
computational  time.  For  the  calculations  shown  here,  we  estimate  that  the  savings  can  be  as  much  as  90%  of  the  total 
computational  time.  Based  on  these  results,  we  note  that  the  ROM  formulation  shows  significant  promise  for  extension 
to  more  realistic  systems  of  equations  such  as  the  Euler  or  Navier-Stokes  equations.  These  extensions  will  be  the 
subject  of  future  studies. 


V.  Conclusions 

A  novel  POD-Galerkin  reduced-order  formulation  is  derived  and  applied  to  a  model  reaction-advection  equation 
that  is  representative  of  combustion  response  problems.  The  approach  utilizes  CFD  solutions  to  obtain  POD  basis 
vectors,  which  are  used  within  a  Galerkin  formulation  to  reduce  the  governing  partial  differential  equations  to  an 
ordinary  differential  equation.  The  resulting  ROM  is  verified  against  known  solutions  and  its  overall  performance 
characteristics  are  investigated.  Specifically,  three  steps  are  taken  in  their  validation:  one,  the  accuracy  of  the  solution 
reconstructed  using  a  certain  number  of  POD  modes  is  assessed;  two,  the  accuracy  of  the  ROM  solution  obtained  by 
solving  the  ODE  is  determined  by  comparing  with  known  CFD  solutions  for  single  and  broadband  forcing;  three,  the 
accuracy  of  the  ROM  solution  to  predict  specific  frequency  response  is  studied  using  the  ROM  that  has  been  trained 
using  broadband  forcing.  These  studies  are  carried  out  for  both  linear  and  non-linear  versions  of  the  model  equation. 
In  addition,  parametric  studies  of  the  effects  of  time-step  size,  discretization  method  and  number  of  POD  modes  on 
the  performance  of  the  ROM  are  also  systematically  evaluated. 

The  results  indicate  that  the  ROM  is  capable  of  predicting  highly  accurate  solutions  with  well-converged  error 
profiles.  It  is  important  that  the  discretization  of  the  ROM  follows  the  spatial  discretization  used  in  the  original  CFD 
solution  procedure.  When  inconsistent  discretization  techniques  are  used,  the  linear  matrix  in  the  Galerkin  formulation 
is  observed  to  have  spurious  eigenvalues  that  prevents  error  convergence.  For  single  frequency  forcing,  two  POD 
modes  are  found  to  be  sufficient  to  capture  nearly  100%  of  the  total  solution,  while  for  broadband  forcing,  a 
concomitantly  larger  number  of  POD  modes  are  needed.  Specifically,  when  nine  discrete  frequencies  were  used  in  the 
forcing,  eighteen  POD  modes  were  required  to  faithfully  represent  the  solution  both  with  the  POD  and  the  ROM, 
which  indicates  that  two  POD  modes  are  needed  for  each  discrete  frequency  involved.  The  ROM  is  also  found  to  be 
effective  for  predicting  discrete  frequency  response  when  the  ROM  has  been  trained  by  broadband  forcing.  This 
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property  is  important  for  practical  application  of  the  ROM  for  design  studies.  It  should  be  noted,  however,  that  the 
predictions  are  poor  when  the  discrete  frequency  is  not  within  the  range  of  the  broadband  frequencies  used  in  the 
training,  which  is  to  be  expected.  Finally,  studies  are  also  carried  out  to  characterize  the  non-linear  response  of  the 
ROM  solution  that  is  trained  using  a  low  (i.e.,  linear)  amplitude  forcing.  Good  predictions  are  obtained  for  low  forcing 
amplitudes,  but  they  become  poor  as  the  forcing  amplitude  becomes  large  enough  for  significant  non-linear  effects. 

Overall,  the  ROM  demonstrates  that  it  can  be  used  effectively  for  carrying  out  parametric  surveys  of  combustion 
response  at  a  fraction  of  the  cost  of  the  full  partial  differential  equation  solution.  Future  studies  will  focus  on  the 
extension  of  the  ROM  approach  to  more  realistic  combustion  physics  involving  the  Euler  and/or  the  Navier-Stokes 
equations. 
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